SIR Modeling Example
The SIR Model
This is a 3-compartment ODE model, in which a person at any one time is either:
- Susceptible to disease
- Infected (assumed infectious)
- Recovered (and thus not infectious)
Image source: https://mc-stan.org/users/documentation/case-studies/boarding_school_case_study.html
Implications
- The total population size is fixed: \(N = S(t) + I(t) + R(t)\) at all times
- The \(R\) state is terminal: once a subject is in this state, the subject’s state cannot change.
- The dynamics of \(I(t)\) over time depend entirely on \(\beta\) and \(\gamma\), specifically:
\[ R_0 = \frac{\beta}{\gamma} \]
This number is called the basic reproduction number, reflecting the expected number of new infections from a single infected person in a population of susceptible people.
src: https://en.wikipedia.org/wiki/Measles
Posterior estimates
Formatting posterior estimates takes some care
Once formatted, we can use tidybayes to summarize flexibly:
Parameters
Data
There is a canonical dataset describing a particularly virulent flu outbreak at a boarding school in 1978. Later, it was determined that the flu was H1N1, however this was not known during the event.
We are given the following in the data description:
The index case was infected by 1978-01-10, and had febrile illness from 1978-01-15 to 1978-01-18. 512 boys out of 763 became ill.
It’s worth noting that the minimal date in our data is: 1978-01-22.
These data are missing details of the index case. We will hard-code this information as our initial conditions.
We will also ignore the “convalescent” category, for now, and lump them in with the “recovered” group (neither infectious nor susceptible).
Posterior estimates
Formatting posterior estimates takes some care
library(tidybayes)
fit_draws <- flu_fit %>%
gather_draws(yrep[i], ys[i, state_index]) %>%
ungroup() %>%
left_join(tibble(t = flu_t, cases = flu_cases) %>% mutate(i = row_number())) %>%
left_join(tibble(state = c('S', 'I', 'R')) %>%
mutate(state_index = row_number())) %>%
mutate(state = if_else(is.na(state), .variable, state)) %>%
select(state, .value, .chain, .iteration, .draw, cases, t) %>%
spread(state, .value)Now, we can compare the posterior estimates to our observed data:
Predictive checks
Parameters
COVID-19 data
Formatting posterior estimates takes some care
library(tidybayes)
fit_draws <- covid_fit %>%
gather_draws(yrep[i], ys[i, state_index]) %>%
ungroup() %>%
left_join(tibble(t = covid_t, cases = covid_cases) %>% mutate(i = row_number())) %>%
left_join(tibble(state = c('S', 'I', 'R')) %>%
mutate(state_index = row_number())) %>%
mutate(state = if_else(is.na(state), .variable, state)) %>%
select(state, .value, .chain, .iteration, .draw, cases, t) %>%
spread(state, .value)
Model formulation
As an example, consider this 8-compartment model used to describe COVID-19 disease dynamics.
The observed data include daily counts of Case, Hospitalization, and Death counts. Each is modeled using a negative binomial given the estimated true prevalence.
src: Keller, JP, Zhou, T, Kaplan, A, Anderson, GB, Zhou, W. Tracking the transmission dynamics of COVID-19 with a time-varying coefficient state-space model. Statistics in Medicine. 2022; 41( 15): 2745- 2767. doi:10.1002/sim.9382
Covariates
In addition, their model includes particular covariate effects on time-varying infectiousness parameters.
For example, one informative covariate was the degree of mobility (from aggregate cell-phone tracking data)
Keller, JP, Zhou, T, Kaplan, A, Anderson, GB, Zhou, W. Tracking the transmission dynamics of COVID-19 with a time-varying coefficient state-space model. Statistics in Medicine. 2022; 41( 15): 2745- 2767. doi:10.1002/sim.9382